#include<Algorithm/mesh_bool.h>

#include<Mesh/libcell_iterator.h>

#define quote lib_cell_quote

// bool运算要用haskell语言的高精度库实现

typedef struct My_Vertex{    
    Mesh* mesh;
    //v1 曲面mesh对应的点
    //v2 相交曲线curve对应的点
    template_v* v1,*v2;
    double* p;
    //store new vertex's id
    int id;
    double dis;
    void* prop;
}My_Vertex;

static inline void my_vertex_init_(My_Vertex* myv)
{
    myv->v1=NULL;myv->v2=NULL;
    myv->mesh=NULL;
    myv->id=-1;
    myv->dis=0;
    myv->p=NULL;
    myv->prop=NULL;
}

static inline  double compute_included_angle(double *a,double*b,double*n)
{
    normalize(a,3);normalize(b,3);
    double*n1=out_product(a,b);
    double re=acos(inner_product(a,b,3));
    if(inner_product(n1,n,3)<0)
    {
        re=2*M_PI-re;
    }
    free(n1);
    return re;  
}


//获取cell的边界点（按顺序）
//tree1是m(c所在mesh)到新建点的映射
//tree2是相交线curve到新建点的映射

//******************
static inline int my_cmp(void *a,void* b)
{   
    return (((My_Vertex*)a)->dis>((My_Vertex*)b)->dis?0:1) ;
}

static  Node*  my_get_cell_boundary_vertices(template_c* c,Mesh2_Crossover_Point*mcp ,
    Int_RB_Tree* tree1,Int_RB_Tree*tree2)
{
    Node* re=NULL;
    // Mesh m;
    // Mesh_init(&m);
    Node* hfs=Mesh_adjust_halffaces(c);
    // Mesh_free(&m);
    for(Node* hfit=hfs;hfit!=NULL;hfit=(Node*)(hfit->Next))
    {
        template_hf*hf=(template_hf*)(hfit->value);
        template_f* f=hf->face;
        template_v* v=hf->vertices[0];
        re=node_pushback(re,tree1->find(tree1,v->id));
        Node* temp_node=NULL;

        for(Node* nit=(Node*)(mcp->f2p->find(mcp->f2p,f->id));nit!=NULL;nit=(Node*)(nit->Next))
        {
            template_v*v1=(template_v*)(nit->value);
            My_Vertex* mv=(My_Vertex*)(tree2->find(tree2,v1->id));
            // if(mv==NULL)
            // {
            //     printf("cuowu\n");
            // }
            mv->dis=distance_of_two_points(mv->v2->point,v->point,3);
            temp_node=node_overlying(temp_node,mv);
        }
        
        
        temp_node= node_bub_sortn(temp_node,my_cmp);
        // if(temp_node!=NULL)
        // {
        //     if(node_size(temp_node)==2){printf("*****************\n");}
        //     printf("begin sort:%d\n",node_size(temp_node));
        //     for(Node* nit=temp_node;nit!=NULL;nit=(Node*)(nit->Next))
        //     {
        //         My_Vertex*mv=(My_Vertex*)(nit->value);
        //         printf("%d dis:%lf\n",mv->v2->id,mv->dis);
        //     }
        // }
        // printf("end sort\n");
        Node* temp_noder=node_reverse(temp_node);
        // printf("end reverse\n");
        re->Next=temp_node;
        if(temp_node!=NULL){temp_node->Prev=re;re=temp_noder;}
    }
    free_node(hfs);

    return node_reverse(re);
}

//返回值在链表的索引
//************
static int my_get_node_value_index(Node* node,void* value)
{
    int re=-1;
    int i=0;
    for(Node*nit=node;nit!=NULL;nit=(Node*)(nit->Next))
    {
        if(nit->value==value)
        {
            re=i;
            break;
        }
        i++;
    }
    return re;
} 
static int array_find_value(void** a,void*v,int len)
{
    for(int i=0;i<len;i++)
    {
        if(a[i]==v)
        {
            return 1;
        }
    }

    return 0;
}



//curve 是相交的曲线
//mcp是m到curve的映射
//vertices储存新建的点，用来释放内存
//给一个cell,返回切割后的多边形区域
static inline int is_same_point(double *p1,double *p2)
{
    if(fabs(p1[0]-p2[0])<1e-6&&fabs(p1[1]-p2[1])<1e-6&&fabs(p1[2]-p2[2])<1e-6)
    {
        return 1;
    }
    return 0; 
}

static void my_insert_one_vertex_faces_to_tree(Int_RB_Tree* tree,template_v*v)
{
    // for(Node* nit=v->faces;nit!=NULL;nit=(Node*)(nit->Next))
    for(int i=0;i<v->facesn.size;i++)
    {
        template_f* f=(template_f*)(v->facesn.extra_data[i]);
        // template_f* f=(template_f*)(nit->value);
        tree->insert(tree,f->id,f);

    }

}

static Node* my_get_split_areas_from_one_cell(template_c*c , Mesh2_Crossover_Point* mcp,
    Mesh* curve,Node** vertices,Node**edges)
{ 
    Int_RB_Tree* tree1=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree1);
    Int_RB_Tree* tree2=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree2);
    Int_RB_Tree*tree_f=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree_f);
    for(int i=0;i<c->vertices_size;i++)
    {
        My_Vertex*mv=(My_Vertex*)malloc(sizeof(My_Vertex));
        my_vertex_init_(mv);
        *vertices=node_overlying(*vertices,mv);
        mv->v1=c->vertices[i];mv->p=c->vertices[i]->point;
        tree1->insert(tree1,mv->v1->id,mv);

        for(Node* nit=(Node*)(mcp->v2p->find(mcp->v2p,c->vertices[i]->id));nit!=NULL;
            nit=(Node*)(nit->Next))
        {
            template_v* v2=(template_v*)(nit->value);
            tree2->insert(tree2,v2->id,mv);
            my_insert_one_vertex_faces_to_tree(tree_f,v2);
        }
      
    }

    //for(Node* hfit=c->halffaces;hfit!=NULL;hfit=(Node*)(hfit->Next))
    for(auto chfit=Mesh_chf_begin(NULL,*c);chfit!=Mesh_chf_end(NULL,*c);chfit++)
    {
        template_f* f=quote(chfit)->face;
        Node* value_n=(Node*)(mcp->f2p->find(mcp->f2p,f->id));
       
        for(Node* nit=value_n;nit!=NULL;nit=(Node*)(nit->Next))
        {
            My_Vertex* mv=(My_Vertex*)malloc(sizeof(My_Vertex));
            my_vertex_init_(mv);
            *vertices=node_overlying(*vertices,mv);
            mv->v2=(template_v*)(nit->value);mv->p=mv->v2->point;

            tree2->insert(tree2,mv->v2->id,mv);
            my_insert_one_vertex_faces_to_tree(tree_f,mv->v2);
        }
    }

    for(Node* nit=(Node*)(mcp->c2p->find(mcp->c2p,c->id)); nit!=NULL;nit=(Node*)(nit->Next))
    {
        My_Vertex* mv=(My_Vertex*)malloc(sizeof(My_Vertex));
        my_vertex_init_(mv);
        *vertices=node_overlying(*vertices,mv);
        mv->v2=(template_v*)(nit->value);mv->p=mv->v2->point;
        // if(c->id==1132)
        // {
        //     printf("id:%d %lf %lf %lf\n",mv->v2->id,mv->v2->point[0],mv->v2->point[1],mv->v2->point[2]  );
        // }
        tree2->insert(tree2,mv->v2->id,mv); 
        my_insert_one_vertex_faces_to_tree(tree_f,mv->v2);
    } 
    
    //printf("begin boundaries\n");
    Node* boundaries=my_get_cell_boundary_vertices(c,mcp,tree1,tree2);
     //printf("end boundaries\n");
    *edges=NULL;
    for(auto it=tree_f->begin(tree_f);it.it!=NULL;it++)
    {
        template_f*f=(template_f*)(it.second);
        My_Vertex**mvs=(My_Vertex**)malloc(sizeof(My_Vertex*)*2);
        mvs[0]=(My_Vertex*)(tree2->find(tree2,f->vertices[0]->id));
        mvs[1]=(My_Vertex*)(tree2->find(tree2,f->vertices[1]->id));
        if(mvs[0]==NULL||mvs[1]==NULL||mvs[0]==mvs[1])
        {
            free(mvs);
        }
        else
        {
            *edges=node_overlying(*edges,mvs);
        }

    }
    int_rb_tree_free(tree_f);
    int_rb_tree_free(tree1);int_rb_tree_free(tree2);
    return boundaries;
} 


//m1是要remesh的网格区域
//curve_inter
//mcp是m1到curve的映射
//tree2 是curve的点到返回mesh的点映射
//tree3 是返回mesh的cell到 m1的cell映射。 

//如果稳定就修改,主要对那些没有交线的cell直接返回。

static inline template_v* my_tool_get_vertex_from_my_vertex(Int_RB_Tree* tree1,Int_RB_Tree* tree2,
    My_Vertex* mv)
{
    return mv->v1!=NULL?((template_v*)tree1->find(tree1,mv->v1->id)):((template_v*)tree2->find(tree2,mv->v2->id));
}

static inline void my_register_new_vertex(template_v*nv,My_Vertex* mv,Int_RB_Tree* tree1,
    Int_RB_Tree*tree2,Mesh2_Crossover_Point*mcp)
{
    if(mv->v1!=NULL)
    {
        tree1->insert(tree1,mv->v1->id,nv);
        for(Node* nit=(Node*)(mcp->v2p->find(mcp->v2p,mv->v1->id));nit!=NULL;
            nit=(Node*)(nit->Next))
        {
            template_v* v2=(template_v*)(nit->value);
            tree2->insert(tree2,v2->id,nv);
        }
    }
    if(mv->v2!=NULL)
    {
        tree2->insert(tree2,mv->v2->id,nv);
    }


}

static void show_subdivision_data(double** v,int len, 
    int* bound,int bound_len,int**edges,int edge_len)
{
    for(int i=0;i<len;i++)
    {
        printf("%lf %lf %lf\n",v[i][0],v[i][1],v[i][2]);
    }
    for(int i=0;i<bound_len;i++)
    {
        printf("%d ",bound[i]);
    }
    printf("\n");
    for(int i=0;i<edge_len;i++)
    {
        printf("%d %d\n",edges[i][0],edges[i][1]);
    }
    printf("end *************\n");
}
//#include "tool_bool.hpp"
Mesh* remesh_along_intersecting_curve(Node* m1,Mesh2_Crossover_Point*mcp,
    Mesh*curve,Int_RB_Tree* tree2,Int_RB_Tree* tree3)
{
    if(m1==NULL)
    {
        return NULL;
    }

    Mesh* re=(Mesh*)malloc(sizeof(Mesh));
    Mesh_init(re);
    re->dimension=2;re->simplex=1;re->manifold_require=1;
    //tree1是m1的点到re的点映射

    Int_RB_Tree* tree1=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree1);
    for(Node* cit=m1;cit!=NULL;cit=(Node*)(cit->Next))
    { 
        Node* vertices=NULL;
        //printf("begin spllit area \n");
        //printf("remsh c id:%d\n",((template_c*)(cit->value))->id);

        Node* n_edges=NULL;
        Node* boundaries=my_get_split_areas_from_one_cell((template_c*)(cit->value),mcp,curve,
            &vertices,&n_edges); 
        int len=0;
        My_Vertex**vertices1=(My_Vertex**)node2array(vertices,&len);
        double**vertices_p=(double**)malloc(sizeof(double*)*len);
        for(int i=0;i<len;i++)
        {
            vertices1[i]->id=i;
            vertices_p[i]=vertices1[i]->p ;
        }
        int bound_len=0;
        My_Vertex**boundaries1=(My_Vertex**)node2array(boundaries,&bound_len);
        int* bound=(int*)malloc(sizeof(int)*bound_len);
        for(int i=0;i<bound_len;i++)
        {
            bound[i]=boundaries1[i]->id;
        }

        int edge_len=0;
        My_Vertex***edges1=(My_Vertex***)node2array(n_edges,&edge_len);
        int **edges=NULL;
        if(edge_len>0)
        {
            //printf("edge len:%d\n",edge_len);
            edges=(int**)malloc(sizeof(int*)*edge_len);
        }

        for(int i=0;i<edge_len;i++)
        {
            edges[i]=(int*)malloc(sizeof(int)*2);
            edges[i][0]=edges1[i][0]->id;
            edges[i][1]=edges1[i][1]->id;
        }

        int trian_len=0;
        // show_subdivision_data(vertices_p,len, 
        //     bound,bound_len,edges,edge_len);

        // printf("begin sub\n");
        // int**tria_ss= subdivision_of_polygon_with_edges1(vertices_p,len, 
        //     bound,bound_len,edges,edge_len,&trian_len);

        // int**tria_ss=optimization_subdivision_of_polygon_with_edges1(vertices_p,len, 
        //     bound,bound_len,edges,edge_len,&trian_len);
        int**tria_ss=optimization_subdivision_of_polygon_with_edges(vertices_p,len, 
            bound,bound_len,edges,edge_len,&trian_len);
        // if(((template_c*)(cit->value))->id==14848)
        // {
        //     assert(false);
        // }
        //printf("end sub\n");
        template_v*vs[3]={NULL};
        for(int i=0;i<trian_len;i++)
        {
            memset(vs,0,sizeof(template_v*)*3);
            for(int j=0;j<3;j++)
            {
                //assert(tria_ss[i][j]<len);
                My_Vertex* mv=vertices1[tria_ss[i][j]];

                vs[j]=my_tool_get_vertex_from_my_vertex(tree1,tree2,mv);
                if(vs[j]==NULL)
                {
                    template_v*nv=re->create_vertexv(re,mv->p,3);
                    my_register_new_vertex(nv,mv,tree1,tree2,mcp);
                    vs[j]=nv;
                }
            }

            if(vs[0]!=NULL&&vs[1]!=NULL&&vs[2]!=NULL)
            {
                template_c*c=re->create_cellv(re,vs,3);
                if(c!=NULL)
                {
                    tree3->insert(tree3,c->id,cit->value);
                }
            }
            free(tria_ss[i]);
        }
        //printf("hrere1\n");
        SAFE_FREE(tria_ss);
        //printf("hrere2\n");
        SAFE_FREE(edges);
        SAFE_FREE(edges1);
        free_node_value(n_edges);
        free_node(n_edges);
        //printf("hrere3n");
        free(bound);
        free(boundaries1);
        free(vertices_p);
        free(vertices1);
        free_node_value(vertices);
        free_node(vertices);
        free_node(boundaries);
    }
    int_rb_tree_free(tree1); 
    return re;
}



//
static Node* get_union_faces_boundaries(Node* node1)
{
    Node* re=NULL;

    while(node1!=NULL)
    {
        Node* node=NULL;
        template_f* f=(template_f*)(node1->value);
        node1=(Node*)(node1->Next);
        node=node_overlying(node,f->vertices[0]);
        node=node_overlying(node,f->vertices[1]);
        int flag=1;
        while(flag==1)
        {
            flag=0;
            template_v* v1=(template_v*)(node->value);
            template_v* v2=(template_v*)(node_reverse(node)->value);
            for(Node* nit=node1;nit!=NULL;nit=(Node*)(nit->Next))
            {
                template_f* f1=(template_f*)(nit->value);
                if(f1->vertices[0]==v1)
                {
                    node=node_overlying(node,f1->vertices[1]);
                    flag=1;node1=node_remove(node1,nit);
                    break;
                }
                else if(f1->vertices[1]==v1)
                {
                    node=node_overlying(node,f1->vertices[0]);
                    flag=1;node1=node_remove(node1,nit);
                    break;
                }
                else if(f1->vertices[0]==v2)
                {
                    node_pushback(node_reverse(node),f1->vertices[1]);
                    
                    flag=1;node1=node_remove(node1,nit);
                    break;
                }
                else if(f1->vertices[1]==v2)
                {
                    node_pushback(node_reverse(node),f1->vertices[0]);
                    
                    flag=1;node1=node_remove(node1,nit);
                    break;
                }
            }
        }
        re=node_overlying(re,node);
    }


    return re;
}

//m要切割的网格
//curve相交曲线
//tree储存curve 到 m 的映射(tree是需要传入,不会修改)(curve存在孤立点，所以tree定义域不完全)

//tree4 储存m新建点到m的映射
static Node* my_adjust_get_moduan_vertices(Mesh* m)
{
    Node* node1=NULL;
    for(auto it=m->v_begin(m);it!=m->v_end(m);it++)
    {
        int size=quote(it)->facesn.size;
        // int size=node_size(quote(it)->faces);
        if(size ==1||size==0)
        {
            node1=node_overlying(node1,quote(it));
        }
    }
    return node1;   
}
//********
// 添加一些faces用来帮正切割的连通性**
//**********
static Node* my_tianjia_threshold_path_faces(Node* node,Int_RB_Tree* tree,Mesh* m,Mesh* curve)
{
    Node* node1=my_adjust_get_moduan_vertices(curve);

    printf("begin my_tianjia_threshold_path_faces %d\n",node_size(node1));
    while(node1!=NULL)
    {

        Node* n1=NULL,*n2=NULL;
        Node* nit=(Node*)(node1->Next);n1=node1;
        template_v* v=(template_v*)(n1->value),*v1=NULL;
        printf("%d %lf %lf %lf\n",v->id,v->point[0],v->point[1],v->point[2]);
        for(;nit!=NULL;nit=(Node*)(nit->Next))
        {
            v1=(template_v*)(nit->value);
            if(distance_of_two_points(v->point,v1->point,3)<1e-4)
            {
                printf("%d %d :dis:%lf\n ",v1->id,v->id,distance_of_two_points(v->point,v1->point,3) ); 
                n2=nit;
                break;
                
            }

        }
        if(n2!=NULL)
        {
            template_v*vv= (template_v*)(tree->find(tree,v->id));
            template_v*vv1= (template_v*)(tree->find(tree,v1->id));
            Node*fvs= get_shortest_path(m, vv,vv1);
            template_v* vs[2]={0};
            for(Node* nit1=fvs;nit1!=NULL;nit1=(Node*)(nit1->Next))
            {
                if(nit1->Next!=NULL)
                {
                    vs[0]=(template_v*)(nit1->value);
                    vs[1]=(template_v*)(( (Node*)(nit1->Next) )->value);
                    template_f*f=m->get_facev(m,vs,2);
                    if(f==NULL){printf("zicuowu\n");}
                    node=node_overlying(node,f );
                }
                //printf("%d\n",((template_v*)(nit->value))->id );
            }
            free_node(fvs);
        }
        node1=node_remove(node1,n1);
        // printf("node size:%d\n",node_size(re));
        node1=node_remove(node1,n2);
    }

    return node;
}


void cut_along_curve_mesh(Mesh* m,Mesh* curve,Int_RB_Tree* tree,
    Int_RB_Tree* tree4,  Int_RB_Tree *tree5)
{
    if(m==NULL)
    {
        return;
    }
    Node* node=NULL,*node1=NULL;

    for(auto fit=curve->f_begin(curve);fit!=curve->f_end(curve);fit++)
    {
        template_v* vs[2]={NULL};
        vs[0]=(template_v*)(tree->find(tree,quote(fit)->vertices[0]->id));
        vs[1]=(template_v*)(tree->find(tree,quote(fit)->vertices[1]->id));
        if(vs[0]==NULL||vs[1]==NULL)
        {
            node1=node_overlying(node1,quote(fit));
            printf("zheli f 修改\n");
            // if(vs[0]==NULL)
            // {
            //     printf("v id:%d\n",quote(fit)->vertices[0]->id);
            // } 
            // if(vs[1]==NULL)
            // {
            //     printf("v id:%d\n",quote(fit)->vertices[1]->id);
            // } 
            continue;
        }

        template_f* f=m->get_facev(m,vs,2);
        if(f==NULL)
        {
            printf("cuowu f vs:%d %d from :%d %d\n",vs[0]->id,vs[1]->id,
                quote(fit)->vertices[0]->id,quote(fit)->vertices[1]->id);
            printf("get_shortest_path\n");
            Node*fvs= get_shortest_path(m, vs[0],vs[1]);
            printf("end get_shortest_path\n");
            for(Node* nit=fvs;nit!=NULL;nit=(Node*)(nit->Next))
            {
                if(nit->Next!=NULL)
                {
                    vs[0]=(template_v*)(nit->value);
                    vs[1]=(template_v*)(( (Node*)(nit->Next) )->value);
                    f=m->get_facev(m,vs,2);
                    if(f==NULL){printf("zicuowu\n");}
                    node=node_overlying(node,f );
                }

                //printf("%d\n",((template_v*)(nit->value))->id );
            }
            free_node(fvs);

        }
        else
        {
            node=node_overlying(node,f);
        }
    }
    printf("begin unios_bounds\n");
    Node* unios_bounds=get_union_faces_boundaries(node1);
    for(Node* nit=unios_bounds;nit!=NULL;nit=(Node*)(nit->Next))
    {   
        Node* node2=(Node*)(nit->value);
        template_v* vs[2]={0};
        vs[0]= (template_v*)(tree->find(tree,((template_v*)(node2->value))->id)) ;
        vs[1]=(template_v*)(tree->find(tree,((template_v*)(node_reverse(node2)->value))->id));
        //printf("unios :%d :%d\n",vs[0]->id,vs[1]->id);
        Node*fvs= get_shortest_path(m, vs[0],vs[1]);
        for(Node* nit1=fvs;nit1!=NULL;nit1=(Node*)(nit1->Next))
        {
            if(nit1->Next!=NULL)
            {
                vs[0]=(template_v*)(nit1->value);
                vs[1]=(template_v*)(( (Node*)(nit1->Next) )->value);
                template_f*f=m->get_facev(m,vs,2);
                if(f!=NULL){node=node_overlying(node,f);}
            }
           // printf("%d\n",((template_v*)(nit->value))->id );
        }
        free_node(fvs);

        free_node(node2);
    }
    free_node(unios_bounds);

    // node=my_tianjia_threshold_path_faces(node,tree, m,curve);
    // //printf("hereend\n");
    // Int_RB_Tree* tree1=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    // int_rb_tree_init(tree1);
    // Int_RB_Tree* tree2=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    // int_rb_tree_init(tree2);
    printf("begin cut along\n");
    Mesh_cut_along_the_curvef(m,node,tree5,tree4);

    free_node(node);
}


